function [] = qmt_preprocess(subjectID, b1_seq, b1_files, t1_seq, t1_files, b0_seq, b0_files, ref_file)
%qmt_preprocess Creates a mask, B1 map, T1 map, and B0 map. 
%   Creates a mask, B1 map, T1 map, and B0 map. 
%
%   Currently supports:
%
%   B1: AFI data
%   T1: VFA data
%   B0: gre_0 data
%
%   **Data files must be in minc format**
%   
%   Author: Mathieu Boudreau
%   Created: March 21st 2013
%

cur_dir = cd;
%% Function input format checks
%

if(strcmp(b1_seq,'afi')&&strcmp(t1_seq,'vfa'))
    
    %Check existence of b1 files
    for ii=1:length(b1_files)
        if(~exist(sprintf('%s/%s_%dd%d_mri.mnc',cur_dir,subjectID, b1_files(ii),ii),'file'))
            beep_sound();
            fprintf('%s_%dd%d_mri.mnc does not exist. Stop.',subjectID, b1_files(ii),ii)
            return
        end
    end

    %Check existence of t1 files
    for ii=1:length(t1_files)
        if(~exist(sprintf('%s/%s_%d_mri.mnc',cur_dir, subjectID, t1_files(ii)),'file'))
            beep_sound();
            fprintf('%s_%d_mri.mnc does not exist. Stop.',subjectID, t1_files(ii))
            return
        end
    end
end
    
if(strcmp(b1_seq,'gre_b0'))
    
    %Check existence of b0 files
    for ii=1:length(b0_files)
        if(~exist(sprintf('%s/%s_%d_mri.mnc',cur_dir, subjectID, b0_files(ii)),'file'))
            beep_sound();
            fprintf('%s_%d_mri.mnc does not exist. Stop.',subjectID, b0_files(ii))
            return
        end
    end   
end


%% Process B1 (AFI) and T1 (VFA)
%


% Open/Create B1 and T1 folder
if(exist(sprintf('%s/b1',cur_dir),'dir'))
    b1_dir = 'b1/';

else
    mkdir('b1')
    b1_dir = 'b1/';
    fprintf('\nb1 folder created\n')
end

if(exist(sprintf('%s/t1',cur_dir),'dir'))
    t1_dir = 't1/';
    
else
    mkdir('t1')
    t1_dir = 't1/';
    fprintf('\nt1 folder created\n')
end

if(strcmp(b1_seq,'afi')&&strcmp(t1_seq,'vfa'))
   ref_file = sprintf('%s_%i_mri.mnc',subjectID, ref_file); 
   fitDataAFI (subjectID, b1_files, t1_files, 't1_map.mnc', 'b1_map.mnc', ref_file)
   system(sprintf('mv b1_map.mnc %sb1_map.mnc',b1_dir))
   system(sprintf('mv t1_map_lin.mnc %st1_map.mnc',t1_dir))
   system(sprintf('mv m0_map_lin.mnc %sm0_map.mnc',t1_dir))

elseif(strcmp(b1_seq,'none')||strcmp(t1_seq,'none'))
    disp('B1 and T1 data not processed.')
    
else
    beep_sound()
    disp('qmt_preprocess can only process AFI and VFA data. Stop.')
    return
end

%% Process B0 (gre_b0)
%

% Open/Create B0 folder
if(exist(sprintf('%s/b0',cur_dir),'dir'))
    b0_dir = 'b0/';
    
else
    mkdir('b0')
    b0_dir = 'b0/';
    fprintf('\nb0 folder created\n')
end

% Call B0 mapping function
if(strcmp(b0_seq,'gre_b0'))
    
    % Split B0 phase images
    if(size(b0_files)==1)
        system(sprintf('mincreshape -clobber  -dimrange time=0,0 %s_%i_mri.mnc  b0_te_1.mnc', subjectID, b0_files))
        system(sprintf('mincreshape -clobber  -dimrange time=1,0 %s_%i_mri.mnc  b0_te_2.mnc', subjectID, b0_files))
    % Copy and rename both B0 phase files if two files were inputed.
    elseif(ischar(b0_files)&&(size(b0_files)==2))
        system(sprintf('cp %s_%s_mri.mnc b0_te_1.mnc',subjectID, b0_files(1)))
        system(sprintf('cp %s_%s_mri.mnc b0_te_2.mnc',subjectID, b0_files(2)))
    else
        beep_sound()
        disp('qmt_preprocess can only process gre_b0 B0 data. Stop.')
        return
    end
    
    calculateB0FieldMap_3t('b0_te_1.mnc', 'b0_te_2.mnc', 'b0_map.mnc')
    system(sprintf('mv b0_map.mnc %sb0_map.mnc',b0_dir))
    
elseif(strcmp(b0_seq,'none'))
    disp('B0 data not processed.')
    
else
    beep_sound()
    disp('qmt_preprocess can only process b0_files of size 1 or 2. Stop.')
    return
end



%% If B0 is 2D and B1/T1 data is 3D, copy B1/T1 files and resample them to match B0;
%

b0_hdr = niak_read_hdr_minc('b0/b0_map.mnc');
b1_hdr = niak_read_hdr_minc('b1/b1_map.mnc');
t1_hdr = niak_read_hdr_minc('t1/t1_map.mnc');

if(or(b0_hdr.info.dimensions(3)<b1_hdr.info.dimensions(3),b0_hdr.info.dimensions(3)<t1_hdr.info.dimensions(3))&&(b1_hdr.info.dimensions(3)==t1_hdr.info.dimensions(3)))
   system('cp mask.mnc mask_3d.mnc')     
   system('mincresample -clobber -like b0/b0_map.mnc -nearest_neighbour mask_3d.mnc mask.mnc') 
   
   system('cp b1/b1_map.mnc b1/b1_map_3d.mnc')
   system('mincresample -clobber -like b0/b0_map.mnc b1/b1_map_3d.mnc b1/b1_map.mnc') 

   
   system('cp t1/t1_map.mnc t1/t1_map_3d.mnc')     
   system('mincresample -clobber -like b0/b0_map.mnc t1/t1_map_3d.mnc t1/t1_map.mnc') 
end

%% Mask files 
% Some masking may have been done during the processing, but this is done
% as a final touch,

system('mincmask -clobber b0/b0_map.mnc mask.mnc b0/temp.mnc')
system('mv b0/temp.mnc b0/b0_map.mnc')

system('mincmask -clobber b1/b1_map.mnc mask.mnc b1/temp.mnc')
system('mv b1/temp.mnc b1/b1_map.mnc')

system('mincmask -clobber t1/t1_map.mnc mask.mnc t1/temp.mnc')
system('mv t1/temp.mnc t1/t1_map.mnc')

 

end


function[] = beep_sound()
    sound([sin(2*(0:0.1:500)),sin(4*(500:0.1:1000)),sin(6*(1000:0.1:1500))]);
end